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Summary 

Nozzle design codes developed earlier under NAG-1-1133 were modified and 
used in order to design a supersonic wind tunnel nozzle with square cross- 
sections. As part of the design process, a computer code was written to 
implement the Hopkins and Hill perturbation solution for the flow in the 
transonic region of axisymmetric nozzles. This technique is used to design 
the bleed slot of quiet-flow nozzles. This new design code is documented in 
this report. 


Development of a Code for Wall Contour Design in the 
Transonic Region of Axisymmetric Nozzles 


by 

Timothy Alcenius (M.S.) and Steven P. Schneider (Assistant Professor) 
School of Aeronautics and Astronautics 
Purdue University 
West Lafayette, IN 47907-1282 


Abstract 

An asymptotic analysis by D.F. Hopkins and D.E. Hill [1] for the flow in the 
transonic portion of an axisymmetric nozzle is used to develop a computer program for 
calculating the nozzle wall contour upstream of the nozzle throat. This code is to be 
integrated with Sivells [2] code for generating axisymmetric nozzles to create a program 
for calculating an entire nozzle contour. The development of the equations for the 
numerical solution is discussed along with a new technique for avoiding singularities. A 
test case was computed, and compared to results from Hopkins and Hill [1] and to results 
from the Douglas Aircraft implementation of Hopkins and Hill [3,4]. 

Introduction 

Many methods are used to design two-dimensional axisymmetric nozzles, one 
being a characteristic method used in a program by Sivells [2]. This code, however, is 
limited when trying to resolve the flow field in the transonic portion of the nozzle. 
Another method was developed by D.F. Hopkins and D.E. Hill [1] in order to address this 
problem. 

Originally, Hopkins and Hill's results were used by Douglas Aircraft Company to 
develop a computer code to calculate the flow in the transonic portion [4]. Even though it 
is still in use today, this code has two problems. First, the code is poorly documented and 
uses logic that is very hard to follow. It is therefore very difficult to understand and use. 
Second, it was programmed when computers were much less powerful than they are now. 
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Even though the solution appears accurate and would run quickly on any machine today, 
this made it even more difficult to understand. 

The purpose of this paper is to describe a new computer program developed to 
implement the results of Hopkins and Hill. The paper includes the development of all the 
relations used by the code and a comparison of the results to the Douglas code and those 
given by Hopkins and Hill. 

The nomenclature used here is identical to that described in Reference 1. 


Development of Code 

In order to use equations (17a-17d) as shown in Hopkins and Hill (also given in 
the appendix). Hr and Mr along with their derivatives with respect to 4 need to be 
calculated. Hopkins and Hill give a relation between Hr and ^ from one dimensional 

considerations along the axis in equation (22) of their paper (see appendix of this paper). 
From this, all of the derivatives were determined (see appendix). Note that 0j is required 
for the determination of M*. Derivatives of Hr(^) up to the fourth order are required for 
03 , along with Mr* derivatives up to second order. These are shown below: 
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where Rs and Ci are parameters that specify the shape of the reference boundary. Since 
the reference boundary and the wall of the nozzle will be of different shapes, Rs and Ci 
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flow field can be calculated for an arbitrarily defined nozzle shape. These functions are 
equations (33) and (34) given in [1] and are repeated in the appendix of this document 

Now that Hr and its derivatives are known, the equations for Mr and its 
derivatives can be obtained. The relations between Mr and Hr are given by equations (13) 
and (16) in Hopkins and Hill (also shown in the appendix). Using the chain rule with 
these equations, the first and second derivatives of Mr can be found. These are given 
below: 



The partial derivatives above are too long to be included here (they are included in the 
appendix). Knowing these variables, any values for ^ and Ti can be used to determine an 
X and y location, flow angle 0, and M*. 

As stated before, any value for q can be chosen to find a solution for equations 
(17a- 17d) by Hopkins and Hill. However, the streamline that ccnresponds to the wall 
contour is the streamline of primary interest Equations (17a-17d) in [1] have been 
normalized so that the value of y at the throat is 1.0. This means that the streamline that 
is desired is the one that passes through this point when the slope of the streamline is 
zero. This occurs when dy/dx = 0, or using the chain rule: 


dx 



= 0. 


(7) 


Since d^dx is not zero in the field. 


( 8 ) 
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is the equation of interest. This equation can be solved explicitly for T| 



with the values for yj and as given in the appendix. 

The equation for ri is then substituted into equation (17b) and y is set equal to 1.0. 
This equation is solved numerically using a bisection method. Then the value for the 
maximum streamline, T), is known from the solution to the above equation as well as the ^ 
location of the minimum point on the wall. This value of ^ is used to shift all of the x 
points so that x=0 will correspond to the minimum point on the wall instead of the 
location of the sonic point along the axis. Although this is not mentioned in Reference 1, 
this procedure is apparently used there also. 

Finally, the initial and final value for ^ must be determined from the user input for 
9a and Xend- The variable 0a represents the maximum angle that the wall makes upstream 
of the throat. The location of this angle is found by iterating over ^ using a bisection 
method to find where the calculated value of 0 is equal to 0a- Since direct computation of 
dy/dx did not agree with (17c) and was consistent with y(x), the calculated value for 9 is 
found from the inverse tangent of dy/dx instead of by using equation (17c). Perhaps this 
mismatch indicates that an accurate description of 0 requires more terms in the series. 
The equation that is used in the program is shown below: 


0 = TAN’’ 




( 10 ) 


where dy/d^ is the same as before and d^dx is given by: 


d^_ 1 

dx 1 - - x'ti'' 


( 11 ) 
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The values for \2 and X 4 ' are given in the appendix. Using the value for rj which 
corresponds to the wall, the value for % which makes 0 equal to 0a can be found. 

The final value for ^ is determined from the input value for Xend- The variable x 
was chosen instead of 0 for use as the ending variable so that a smooth transition to 
Sivells code [2] could be obtained. The value for ^ that corresponds to Xgnd at the wall 
for any other streamline is found using a bisection iteration scheme. 

The equations must be solved at each point in the streamwise direction. Instead of 
using the streamwise coordinate ^ as the independent variable, Mr* is used. This is so the 
iteration to find Mr* only needs to occur to find the beginning and ending locations in the 
calculation. Given Mr* at each location, all of the other variables can be solved for 
explicitly. 

Difficulties were experienced due to singularities in M/ and Mr" at the sonic 
point. These were due to the fact that dMr*/dHr goes to infinity and Hr’ goes to zero at 
Mr*=l (see equation (5)). In order to avoid this problem, a penurbation solution was 
used in the vicinity of Mr*=l. This solution gives: 


dH, ' V(y+i)Rs 


( 12 ) 


where R$ is the same as before and is shown in the appendix. This equation yields results 
that are smooth and correct in the region of Mr* = 1. A perturbation solution was also 
used for Mr '. This was derived by rewriting the equation as: 
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Since in the area of Mr* = 1, dMr*/d^ is constant, it is assumed that the derivative with 
respect to ^ of this function is zero. This leaves: 
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( 14 ) 



which yields a continuous variation of Mr " in this region. 

The final decision that had to be made was the region in which to apply these 
solutions. As can be seen in figure #1, Mr* = 1±.01 gives reasonable results. 


Results 

The validity of the results were checked in two different ways. First, the plots that 
are given by Hopkins and Hill were reproduced. These can be seen in figures 2 to 4. The 
data appear correct, but a more precise comparison was desirable. 

In order to make a second check, results from the Douglas program [3,4] were 
used. Given the variables shown below, the streamline that corresponds to the wall 
contour were compared. 

Table #1 - Parameters Used in Comparison with Langlev Results 


Ra 


0a 

450 

Y 

1.4 


The results are shown in figures 5 to 8. As can be seen in these figures, the results 
from the new code appear to match precisely with the results given by the Douglas code 
[3,4]. By checking the numerical results from both methods, it was found that the codes 
always agree to at least four figures for all variables except for 6 which agrees to two 
figures. Figure 8 also shows that the solution given by the new program is more 
continuous near M* = 1 than the Douglas code [3,4]. Therefore, the new program yields 
a slightly more accurate solution than was previously available. 
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Conclusions 

A new FORTRAN-77 computer code has successfully been developed for 
calculating the transonic flow field in an axisymmetric nozzle. The results of this code 
are in excellent agreement with results from a previous FORTRAN-4 code. The code is 
well documented and commented to allow ready use. The algorithm for the code is also 
provided so that the user may have a better understanding of the equations and the logic 
behind the code development. 

This code has also been integrated with the Sivells [2] code. The integrated 
program is capable of generating an entire nozzle contour, including a region upstream of 
the throat. Although Sivells uses Hall's transonic solution downstream of the throat [2], 
the two asymptotic solutions should agree very near the throat. Current results, to be 
reported on elsewhere, indicate that this is the case. 
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Appendix 

The following equations were taken directly from Hopkins and Hill (1). They are 
placed here fore the readers convenience. 
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1(0.1277-0,0355(Y'-1.4)1 


arctan 


lnl5000-lnR, 
19.55 + 1.25(7-1.4) 


Here, Ra and 0a are the actual values for the radius of curvature and maximum inflow 
angle that are input by the user. 

The following equations were derived for use in the program. 

e; = 1 {h;"(h,h;mJ + h^m^m;) + h;{(h;)^mJ + h^h;m^ +4h,h;m,m; + 
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Figure #2- Re-computation of Figure 4 from Hopkins and Hill [1] 
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Figure #5-Comparison of X between Current Program and Douglas Program 
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Appendix 

An electronic version of the source for the Fortran-77 code can be obtained 
from the Experimental Methods Branch at the NASA Langley Research 
Center, Hampton, VA, 23681. Contact Dr. Stephen Wilkinson. Alter- 
natively, contact Prof. Steven P. Schneider, Aeronautical and Astronauti- 
cal Engineering, Purdue University, West Lafayette, IN 47907-1282, email 
stevesSecn.purdue.edu, telephone (317) 494-3343. 

The following pages contain a printed version of the current source code. 



* This is a program to test the hopkins-hill subroutine sps 6-2-93 

* subroutine hophill (ra, thetaa,gam,thetae,yt,n,xs,ys,amachs) 
real xs (100) ,ys (100) , amachs (100) 

data n/100/ 

* 

print *, 'What is the radius of curvature?' 
read *, ra 

print *, 'What is the initial angle (degrees)?' 
read * , thetaa 

print *, 'What is the ratio of specific heats?' 
read * , gam 

print *, 'What is the final value of x?' 
read *, xexit 

* 

write (*,*) 'opening hophill.tst for test output' 
open (unit=l, file='hophill . tst' ) 
write (1,*) 'x,y,amach along contour' 
y t = 1.0 

call hophil 1 ( ra , thetaa , gam , xexit , yt , n , xs , ys , amachs) 
write (*,*) 'returned from first call to hophill' 
do 50 i = l,n 

write ( 1 , * ) xs ( i ) , y s ( i ) , amachs ( i ) 

50 continue 
yt = 0.95 

write (1,*) 'x,y,amach along streamline through yt= ',yt 
cal 1 hophill ( ra , thetaa , gam , xexit , y t , n , xs , ys , amachs ) 
do 80 i = l,n 

write (1,*) xs(i) ,ys(i) , amachs (i) 

80 continue 

stop 'normal end' 
end 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


From alcenius Wed Jun 2 02:58:29 1993 

This program is developed to calculate the transonic flow 
region of an axisymmetric nozzle. The equations are taken from 
Hopkins, D.F. and Hill, D.E., "Effect of Small Radius of 
Curvature on Transonic Flow in Axisymmetric Nozzles", AIAA 
Journal Vol.4, PP. 1337-1343 (August 1966). 

Written by Tim Alcenius 5-24-93. 

The variables in the main program are defined as follows: 

ra - The actual radius of curvature input by the user 
thetaa - The actual initial flow angle of the wall input by 
the user. 

xend - The last value of x that data is written out for, 
input by the user. 

gam - The ratio of specific heats input by the user, 
cl - Reference bounday shape parameter for flow angle Eqn. 34. 
rs - Reference bounday shape parameter for radius of curvature 
Eqn . 33. 

cltop - The numerator of Eqn #34 
clbottom - The denominator of Eqn #34 
gmp - Gamma + 1 
gmm - Gamma - 1 


* The first part of the program defines all of the variables and 

* reads in the user directed input for the problem. 

* This is a subroutine version of the same program, modified to 

* allow finding interior streamlines, sps 6-2-93 


subroutine hophi 1 1 ( r a , thetaa 1 , garni , xendl , yt 1 , n , xs , ys , amachs ) 
parameter (ndim=2) 

* xs,ys are arrays with coords of line, amachs is mach no. at pt, n is no. pts. 

* yt is the throat value of y for the streamline to be returned 

real ra, thetaa, gam, gmm, gmp, cltop, clbottom, cl, rs, mst 
real xend, pi, mrst, mrstthroat, mrstbeg,mrstend,rorstnow,mrststr 
real xs(n) ,ys(n) , amachs (n) 

real stpts (ndim+1 , ndim) ! initial simplex guess for AMOEBA 
real stans (ndim+1) lvalues of GETSTR at three pts in stpts 
real guess(ndim) 1 for guesses 

logical linit , Isuperl 1 initialization variable, first call, yt=l 
save linit 

external amoeba, getstr, getxy 

common/stream/yt Ito pass the desired value of yt to getstr 
common/ input/ rs, cl, gam 1 input values for passing to subroutines 
common/ throat/etaw,xithroat, mrstthroat, xslO ! results at wall, throat from 

* first iteration for continuation. Note that x,y = 0,1 at throat! 

* remember that mrst depends only on xi, not on eta! 

common/ends/mr stbeg , mrstend , thetaa , xend , thetanow , mrstnow , xbeg 

* Imrst at entry and exit, set in MRSTBAE, theta at entry and exit 

data tiny/1. Oe-5/ ! allowable error in iteration 

data ftol/l.Oe-6/ ! allowable error in AMOEBA iteration 


thetaa = thetaa 1 
xend = xendl 
yt = ytl 

gam = garni !to avoid fortran error, formal argument illegal in common 
if (yt .eq. 1.0) then ! initialize 
if (linit) then 

write (*,*) 'should only call with yt=l the first time!' 
stop 'fatal error' 



else 

linit = .true, 
end if 

pi = AC0S(-1.) 

* write (*,*) 'What is the radius of curvature?' 

* read *, ra 

* write (*, *) 'What is the initial angle (degrees)?' 

* read *, thetaa 

* write (*,*) 'What is the ratio of specific heats?' 

* read *, gain 

gnun = gam - 1. 
gmp = gam + 1. 

* This part of the program calculates the top and bottom portions 

* of egn #34 and the calculates the reference boundary shape 

* parameters using Eqn #33 and #34 from the paper. 

cltop = 180. *ATAN( (LOG (thetaa) /(LOG ( 15000.0) -LOG (ra) )) )/pi 
clbottom = 19.55 + 1.25* (gam - 1.4) 

cl = (cltop/clbottom) **(!./ (0.1277 - (gam-1.4) *0. 0335) ) 
rs - (0.9322 + 0. 0565*gam) * (0. 6173*EXP(0. 48/cl** (0. 189) ) + 

+ ra*(1.1234 + 0.00771*cl - 0 . 000163*cl**2) - ra**2* (0 . 0182 + 

+ 0.00111*cl - 0.0000201*cl**2) ) 

* This part sets the angles to radians and calls the program to 

* calculate the flow field properties in the throat area. 

thetaa = -thetaa*pi/180 !here correct signs so input>0 
call mstar Isets up mrstbeg and mrstend 

* write (*,*) 'HOPHILL: calling getxy with throat conditions' 

* write (*,*) 'mrstthroat= ' ,mrstthroat, ' etaw= ',etaw 

call getxy(mrstthroat,etaw,xi,x,y,mst) Ireturns xi, x, y, and mst 

* for guess of mrst and eta 
if (abs (xi-xi throat) .gt. tiny) then 

write (*,*) 'HOPHILL: xi error in getxy' 
end if 

if (abs(x).gt. tiny) write (*,*) 'x error in getxy' 
if (abs(y-l. 0) .gt. tiny) write (*,*) 'y error in getxy' 
write (*,*) 'HOPHILL: proceeding to compute contour pts' 

* Now have the setup for the streamline through the y=l pt. Now 

* get the coords along the streamline, use the subroutine 

delmrst = (mrstend-mrstbeg) /float (n-1) 
do 100 i = l,n 

mrst = mrstbeg+delmrst* ( i-1) 

call getxy (mrst, etaw,xi,x,y, mst) Ireturns xi,x,y,mst for mrst and eta 
xs(i) = X 
ys(i) = y 

amach = sqrt( 2.0*mst**2/( gmp - gmm*mst**2) ) !eqn 16 
amachs(i) = amach 
100 continue 

return !from first call 

else Icall for streamline, use same xbeg and xend as before! 

* Note that lines of constant eta are streamlines - see eqn (9) in paper 

if (.not. linit) then 

write (*,*) 'should call with yt=l the first time!!' 
stop 'fatal error' 
end if 

* Here, iterate for the values of eta and xi at the throat 

* for the streamline selected 

* Do this using minimization of (x**2 + (y-yt)**2), Numerical Recipes 


* routine amoeba, see also Hildebrand. Function GETSTR implements part. 

* Need good guesses. Note that eta approx. = y near throat. Basis 

* of second guess. 

stpts(l,l) = mrstthroat 

stpts(l,2) = etaw ! throat is first pt, space is (mrst,eta) 
stpts(2,l) = 1.0 ! sonic pt near x=0,y=yt 

stpts(2,2) = yt 

stpts(3,l) = 0.95 !just upstream 
stpts(3,2) = yt 
do 280 i = l,ndim+l 
guess (1) = stpts(i,l) 
guess (2) = stpts(i,2) 
stans(i) = getstr (guess) 

280 continue 

call amoeba (stpts, stans, ndim+1, ndim,ndim, ftol , getstr, iter) 
write {*,*) 'HOPHILL: AMOEBA returns after ',iter,' iterations' 
write (*,*) 'values at three points are: (yt= ',yt,')' 

write (*,*) ' (mrst,eta,x,y,abs(x)+abs(y-yt) ,xi) ' 
toterror = 0.0 
do 290 i = 1, ndim+1 

call getxy( stpts (i, 1) , stpts (i, 2) ,xi,x,y,mst) 

* guess (1) = stpts (i,l) 

* guess (2) = stpts (i, 2) 

* zeroit = getstr (guess) 

* write (*,*) 'HOPHILL: guess, zeroit= ', guess, zeroit 
write (*,288) stpts(i, 1) ,stpts(i,2) ,x,y,stans(i) ,xi 

288 format(lx,6(fl2.7,lx) ) 

toterror = toterror+stans (i) 

290 continue 

if (toterror .gt. tiny) then 

write (*,*) 'HOPHILL: toterror is ', toterror 
write (*,*) 'converged to poor solution from AMOEBA' 
stop 'not worth continuing' 
end if 

mrststr = (stpts (1 , 1) +stpts (2 , 1) +stpts (3 , 1) )/3 . 0 
etastr = (stpts (1, 2) +stpts (2 , 2) +stpts (3 , 2) )/3 . 0 

* Now get New mrstbeg and mrstend for this new streamline, using xbeg and xend 

Isuperl = .false. 

mrstbeg = findmrst (etastr , xbeg, Isuperl) ! function for this purpose 
Isuperl = .true. !mrst supersonic at throat or downstream (Not mst!) 
mrstend = findmrst (etastr , xend, Isuperl) 

* 

delmrst = (mrstend-mrstbeg) /float (n-1) 

write (*,*) 'mrststr , etastr , delmrst= ', mrststr, etastr, delmrst 
do 300 i = l,n 

mrst = mrstbeg+delmrst* (i-1) 

call getxy (mrst, etastr, xi,x,y, mst) Ireturns xi,x,y,mst for mrst, etastr 
xs(i) = X 
ys(i) = y 

amach = sqrt( 2.0*mst**2/( gmp - gmm*mst**2) ) !eqn 16 
amachs(i) = amach 
300 continue 
return 
end if 
end 

* FUNCTION GETSTR 

* this function serves to interface between subroutine GETXY and 

* the Numerical recipes routine AMOEBA 

* 


function getstr (guess) 

real guess (2) ,mrst ,mst ! (mrst and etast) 

coiranon/stream/yt ! desired streamline is through x=0,y=yt 

* 

mrst = guess (1) 
etast = guess (2) 

call getxy (mrst , etast , xi , x , y , mst ) 

getstr = abs(x) + abs(y-yt) Iminimize this, will have soln 

* write (*,1) x,y, getstr, yt 

* 1 format (' GETSTR: x,y , getstr, yt= ' , 4 (fl3 .7, lx) ) 

* write (*,*) 'mrst,etast= mrst, etast 

return 

end 

* FUNCTION GETX2 

function getx2(mrst) 
real rorst,mstout 

* given mrst, passed values, calls getxy to get x, compares to desired value 

common/getx2pass/etapass , xpass 

* 

call getxy (mrst , etapass , xiout , xout , yout , mstout ) 

getx2 = xout-xpass 

return 

end 

* FUNCTION FINDMRST 

* finds the mrst on a given streamline given x 

function f indmrst (etastr , x, Isuperl) 

real mrst ,mrstlow,mrsthigh,mrstmax,mrstbeg,mrstend 
logical succes, Isuper, Isuperl 
external getx2 

common/ends/mrstbeg , mr stend , thetaa , xend , thetanow , mr stnow , xbeg 

* Imrst at entry and exit, set in MRSTBAE, theta at entry and exit 
common/ sonic/ Isuper 

common/getx2pass/etapass , xpass 

common/ input/rs , cl , gam ! input values for passing to subroutines 
data eps/1 . Oe-5/ , small/0 . 3/ 

* 

mrstmax = sqrt ( (gam+1. 0)/ (gam-1. 0) ) 

Isuper = Isuperl 

etapass = etastr ! passing to getx2 
xpass = X 
if (Isuper) then 
mrstlow = 1.0+eps 

mrsthigh = mrstend+small Ismail change from last value, 

* if value too big, problem with h's 

if (mrsthigh .gt. mrstmax-eps) then 

write (*,*) 'FINDMRST: mrstend is high- ',mrstend 
mrsthigh = mrstmax-eps 
end if 

else . . 

mrstlow = mrstbeg I else will have problem with h's in itermrst 

mrsthigh = 1.0-eps 
end if 

xzlow = getx2 (mrstlow) 
xzhigh = getx2 (mrsthigh) 
if (xzlow*xzhigh .It. 0.0) then 
succes = .true, 
else 

succes = 


. false. 


end if 

* remove zbrac, goes off too far 

* call zbrac(getx2,mrstlow,mrsthigh,succes) 
if (succes) then 

write (*,*) 'FINDMRST: bracketed x with mrst in ' , 
> mrstlow,mrsthigh 

else 

stop 'FINDMRST: failed to bracket x' 
end if 

mrst = rtbis(getx2,mrstlow,mrsthigh,eps) 

* make one more call to set up final params 

xzero = getx2(mrst) Hast call to check 
write (*,*) 'FINDMRST: x conv. within: ', xzero 
findmrst = mrst 

* 

return 

end 


* SUBROUTINE MSTAR: 
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The purpose of this subroutine is to calculate Eqn #17 given 
a specific value of mrstar on the axis and a value for a 
streamline eta. The maximum value of eta is found by calling 
the subroutine etawall. If igetgrid is 1, then 
makes an array of data for contour plots: 

Then, at each point xi along the axis, 

the values for hr and mr and their derivatives are calculated. 

Using this information, eta is incremented from the axis (eta=0) 
to the wall value (eta=etaw) . At each one of these points, the 
corresponding value for x, y, theta, and mrstar is calculated. 

These are the variables used in this subroutine; 

hr - The value of h along the reference line (nozzle axis) . 

Defined in Eqn. #12. 

mr - The value of the Mach nuber along the reference line. 

&&pr - The primed value of the function (either hr, mr, mrst, or thetal) 

&&dpr - The second derivative of the function (either hr or mr) 

hrtpr - The third derivative of hr 

hrfpr - The fourth derivative of hr 

xi - Velocity potential with units of length 

etaw - The value of eta which corresponds to the wall 

eterm - Epxonential of Eqn #22 

x2 - Second term in series for x. Defined in Eqn #lla, #17a 

x4 - Fourth term in series for x. Defined in Eqn Hla, #17a 

yl - First term in series for y. Defined in Eqn #llb, #17b 
y3 - Third term in series for y. Defined in Eqn #llb, #17b 
thetal - First term in series for theta. Defined in Eqn #llc, #17c 
thetas - Third term in series for theta. Defined in Eqn #llc, #17c 
q2 - Second term in series for mstar. Defined in Eqn #lld, #17d 

q4 - Fourth term in series for mstar. Defined in Eqn #lld, #17d 

eta - The value of eta used from eta ~ 0 to eta = etaw. 
mrst - The value of mstar on the reference line, 
ra - The actual radius of curvature input by the user 
cl - Reference bounday constant as defined in Eqn #34 
rs - Reference bounday radius of curvature as defined in Eqn #33 
thetaa - The actual initial flow angle of the wall input by 
the user. 

dmrdmrst - Derivative of mr wrt mrstar for mrpr calculation 
dmrstdhr - Derivative of mrstar wrt hr for mrpr calculation 



* mrstbeg - The initial value of mrst determined from thetaa in mrstbae 

* mrstend - The final value of mrst determined from xend in mrstbae 

* dmrst - Delta mrst for the iteration in the direction of constant xi 

* eden - The denominator for the exponential term 

* deleta - Delta eta for the iteration in the direction of constant eta 

* ddmrdmrst - Derivative of dmrdmrst wrt xi for mrdpr calculation 

* ddmrstdhr - Derivative of dmrstdhr wrt xi for mrdpr calculation 

* xislO - The value of xi where dy/dx = 0 for shifting plots so minimum 

* point on top streamline is located at xi = 0. 

* ddcl - Constant for calculation of mr derivatives 

* ddc& - Constants 2-4 for calculation of ddmrstdhr 

* thc& - Constants 1-5 for calculation of thetalpr 

* imax - The maximum number of points in the xi direction 

* jroax - The maximum number of points in the eta direction 

* gam - The ratio of specific heats 

* gmm - gam - 1 . 

* gmp - gam + 1 . 

* 

* Define the variables used in this subroutine and call etawall 

* in order to find the eta which corresponds to the wall value. 

* 

subroutine mstar 

parameter (igetgrid=0, imax = 2, jroax = 2) 

* only fill grid if igetgrid=l, then make imax,jmax nonzero 

real cl, rs, gmm, gmp, gam 

real mrstbeg, mrstend, mdc4pr, dmrst, mrstnow 

real xi, hr, hrpr, hrdpr, hrtpr, hrfpr, deleta, q4 , etaw 

real mr, dmrdmrst, xend, thetaa, mrst, eta, eterm, eden 

real ddcl, dmrstdhr, mrpr, mrstpr, ddmrdmrst, mdc2, rodc3 

real mdc4, mdcl, mdclpr, mdc2pr, mdc3pr, mrdpr, thcl, thc2 , thc3 

real thc4, thc5, theta3pr, thetal, theta3, x2, x4, yl, y3, q2 

real x(imax,jmax) ,y (imax, jroax) , theta (imax, jmax) ,machst ( imax, jmax) 

real thclpr, thc2pr, thc3pr, thc4pr, thcSpr, thc6pr 

real xslO, label (8), mrstthroat , mrstmax 

integer i, j 

common/input/rs , cl , gam ! input values for passing to subroutines 
common/throat/etaw, xithroat, mrstthroat, xslO iresults at wall, throat from 

* first iteration for continuation. Note that x,y = 0,1 at throat! 

* remember that mrst depends only on xi, not on etal 

common/ ends/mrstbeg , mrstend , thetaa , xend , thetanow , mrstnow , xbeg 

* 

gmm = gam-1.0 
gmp = gam+1.0 
if (igetgrid .eq. 1) then 
open(2, f ile='hophill . in' ) 

write(2,*) 'TITLE=”Countours in Throat Region"' 
write(2 , *) 'VARIABLES="X" , "Y" , "THETA" , "M*" , "XI" , "ETA" ' 
write (2 , *) 'ZONE T=" Contours" , 1=' , imax, ' , J=' , jmax, ' , F=POINT' 
end if 

call etawall ! returns values in the common block param (etaw, xslO) 

write (*,*) 'Eta for the wall =',etaw 

call mrstbae ! returns values in common block /ends/ 

write (*,*) 'The initial Mr* =', mrstbeg 

write (*,*) 'The final Mr* =', mrstend 

if (igetgrid ,eq. 1) then 

dmrst = (mrstend - mr s tbeg )/ FLOAT ( imax- 1) 
if (dmrst .le. 0.0) stop 'MSTAR: dmrst eq 0, fatal' 

* 

* Set the number of points in the eta direction to be 100. This can 


be changed by changing the denomenator on the deleta term. 

Then set constants and the initial integer counter in the x direction. 
Finally, enter the do loop to move along in the x direction on 
axis (eta = 0) . 

deleta = etaw/FLOAT ( jmax-1) 
eden = 2.*rs*cl 
z = . 5 

do 5 k = 1,8 
label (k) = z 
z=z+. 1 
continue 

mrstmax = sqrt (gmp/gmm) 
do 10 i = l,imax 

mrst = mrstbeg+ (i-1) *dmrst 

if (mrst .le. 0) stop 'MSTAR: mrst le 0, fatal' 

if (mrst .ge. mrstmax) stop 'MSTAR; mrst ge mrstmax, fatal' 

Calculate hr and its derivatives and mr and its derivatives. Also 
calculate the value of xi corresponding to the given value for 
mrstar from the do loop. 

hr = SQRT ( (1. /mrst) * (0. 5*gmp - 0. 5*gmro*mrst**2) ** (-1./ (gmm) ) ) 
xi = SQRT(2.*rs*cl*LOG(cl/(cl - hr + 1.))) 
if (mrst .LE. 1.) then 
xi = -xi 
endif 

eterm = EXP(-xi**2/eden) 
hrpr = (xi/rs) *eterm 

hrdpr = (1 ./rs) *eterm - (xi**2/ (rs**2*cl) ) *eterro 

hrtpr = -(3 . *xi/ (rs**2*cl) ) *eterm + (xi**3 ./ (rs**3*cl**2) ) * 

+ eterm 

hrfpr = -(3 ./ (rs**2*cl) ) *eterm + (6. *xi**2/ (rs**3*cl**2) ) * 

+ eterm - (xi**4/ (rs**4*cl**3) ) *eterm 

mr = SQRT(2*mrst**2/ (gmp - gmm*mrst**2) ) 
ddcl = 0.5*gmp - 0. 5*gmm*mrst**2 

dmrdmrst = (2 . *mrst*gmp/ (mr* (gmp - gmm*mrst**2) **2) ) 
if (mrst .LT. 0.99 .OR. mrst .GT. 1.01) then 
dmrstdhr = 2 ./ (mrst*hr/ddcl - hr/mrst) 
mrpr = dmrdmrst*dmrstdhr*hrpr 
mrstpr = mrpr/ dmrdmrst 

ddmrdmrst = (4 . *mrstpr*gmp*ddcl - 4 . *mrst*gmp* (mrpr*ddcl - 
+ 2 . *mr*gmm*mrst*mrstpr) ) / (mr**2* (2 . *ddcl) **3) 

mdcl = gmp*mr 
mdclpr = gmp*mrpr 
mdc2= mrst**2 
mdc2pr = 2 . *mrst*mrstpr 
mdc3 = hr 
mdc3pr = hrpr 
mdc4 = hr/mr**2 

mdc4pr = (hrpr*mr - 2 . *hr*mrpr)/ (mr**3) 

mrdpr = ( (mdclpr* (mdc2* (mdc3-mdc4) ) -mdcl* (mdc2pr* (mdc3-mdc4 ) + 

+ mdc2* (mdc3pr-mdc4pr) ) )/ ( (mdc2* (mdc3-mdc4) ) **2) ) *hrpr+ 

+ dmrdmrst *dmrstdhr*hrdpr 

else 

mrpr = dmrdmrst*SQRT (2 ./ (gmp*rs) ) 
mrstpr = mrpr/ dmrdmrst 

ddmrdmrst = (4 . *mrstpr*gmp*ddcl - 4 . *mrst*gmp* (mrpr*ddcl - 
+ 2 . *mr*gmm*mrst*mrstpr) ) / (mr**2* (2 . *ddcl) **3) 

mrdpr = ddmrdmrst*SQRT(2 ./ (gmp*rs) ) 


endif 


Calculate the values for thetal, theta3, thetaSpr, yl, y3, x2 , 
x4, q2, q4 from Eqn #17a-#17d. All of these remain constant at 
any eta for a given value of xi since they are only dependant 
on hr, mr, and their derivatives. 

thcl = hrdpr 
thclpr = hrtpr 
thc2 = hr*hrpr*mr**2 

thc2pr = (hrpr*mr)**2 + hr*hrdpr*mr**2 + 2*hr*hrpr*mr*mrpr 
thc3 = hr**2*mr*mrpr 

thc3pr = 2*hr*hrpr*mr*mrpr + (hr*mrpr)**2 + hr**2*mr*mrdpr 
thc4 = hr**2*hrtpr 

thc4pr = 2*hr*hrpr*hrtpr + hr**2*hrfpr 
thc5 = mr**2-l 
thcSpr = 2*mr*mrpr 
thc6pr = hrpr**2*hrdpr 

theta3pr = 0 . 25* (thclpr* (thc2+thc3) + thcl* (thc2pr+thc3pr) ) + 
+ 0 . 125* (thc4pr*thc5 + thc5pr*thc4 + thcSpr) 

thetal = hrpr 

theta3 = (1 ./4 .)* (hrdpr* (hr*hrpr*mr**2 + hr**2*mr*mrpr) ) + 

+ (1./8. ) *hr**2*hrtpr*(mr**2-l. ) + (1./24 . ) *hrpr**3 

x2 = 0.5*hr*hrpr 

x4 = (3 ./32) * (hr**2*hrpr*hrdpr* (mr**2-l) ) - (1 ./96 . ) * (hr* 

+ hrpr**3) + theta3*hr*0 . 25 

yl = hr 

y3 = (hr/8) * (hr*hrdpr* (mr**2-l) - hrpr**2) 
q2 = 0.5*hr*hrdpr 

q4 = 0 . 25* (hr**2*hrdpr**2 ) + ( 3 ./32 . ) * (hr**2*hrdpr**2* (mr**2- 
+ 1) ) + (1/32 . ) * (hr*hrdpr*hrpr**2) + 0 . 25*hr*theta3pr 

x2pr = (0.5* (hrpr**2 + hr*hrdpr) ) 

x4cl = (3 ./32 . ) * ( (mr**2-l) * (2 . *hr*hrpr**2*hrdpr +hr**2* 

+ hrdpr**2+hr**2*hrpr*hrtpr) +2 . *hr**2*hrpr*hrdpr*mr*mrpr) 

x4c2 = (1 ./96 . ) * (hrpr**4 + 3 . *hr*hrpr**2*hrdpr) 
x4c3 = 0. 25* (theta3pr*hr + theta3*hrpr) 
x4pr = (x4cl - x4c2 + x4c3) 
ylpr = hrpr 

y3pr = (hrpr/8 . ) * (hr*hrdpr* (mr**2-l) - hrpr**2) + (hr/8.)* 

+ (hrpr*hrdpr* (mr**2-l) + hr*hrtpr* (mr**2-l) + 2*hr* 

+ hrdpr *mr*mrpr - 2*hrpr*hrdpr) 


Set the counter value for the arrays for different eta (0->etaw) . 

Then enter the do loop for the calculation of Eqn #17. 

eta = 0. 

do 20 j = 1, jmax 

Calculate Eqn #17a. Then, if mrstar is less than one, make the value 
for X the negative value for plotting purposes. This is because 
the value for xi is always positive along the axis. However, since 
the paper makes the assumption that the sonic line on the axis occurs 
at X = 0 (pp. 1339 paragraph before section #3, and figure #3) , it 
is known that for all values of mrstar less than 1, the x value 
should be negative. 

x(i,j) == xi - x2*eta**2 - x4*eta**4 

Shift the X axis by xslO so that the point where x = 0 corresponds 


* 

* 


to the nozzle throat. 


x(i,j) = x(i,j) - xslO 

* 

* Calculate Eqn #17b. 

* 

y(i,j) = yl*eta + y3*eta**3 

* 

* Calculate Eqn #17c. 

* 

theta (i,j) = ATAN ( (ylpr*eta+y3pr*eta**3) / (1 . -x2pr*eta**2 - 
+ x4pr*eta**4) ) 

* 

* Calculate Eqn #17d. 

* 

inachst(i,j) = mrst*(l + q2*eta**2 + q4*eta**4) 

* 

* Change counters for the next values. 

* 

write(2,99) x(i, j) ,y(i, j) ,theta(i, j) , machst ( i , j ) ,xi,eta 
eta = eta + deleta 
20 continue 

10 continue 

close (2) 
end if 

99 Format(6(lx,el4.7) ) 

return 
end 


* SUBROUTINE ETAWALL (name changed from etamax sps 6-2-93) 

* This subroutine is used to find the streamline which will 

* correspond to the contour for the wall of the nozzle. This is 

* done using a bisection method to iterate of xi to find the 

* value of eta where dydx = 0 and y = 1 as defined by the paper 

* in section 4, pp 1340. 

* 

* The variables that are used in this subroutine are listed below: 

* 

* hr - The value of h along the reference line (nozzle axis) . 

* Defined in Eqn. #12. 

* mr - The value of the Mach nuber along the reference line. 

* &&pr - The primed value of the function (either hr, mr, yl, or y3) 

* hrdpr - The second derivative of hr 

* hrtpr - The third derivative of hr 

* xi - Velocity potential with units of length 

* etaw - The value of eta which corresponds to the wall 

* eterm - Epxonential of Eqn #22 

* eden - The denominator for the exponential term 

* x2 - Second term in series for x. Defined in Eqn #lla, #17a 

* x4 - Fourth term in series for x. Defined in Eqn #lla, #17a 

* yl “ First term in series for y. Defined in Eqn #llb, #17b 

* y3 - Third term in series for y. Defined in Eqn #llb, #17b 

* theta3 - Third term in series for theta. Defined in Eqn #llc, #17 

* mrst - The value of mstar on the reference line. 

* cl - Reference bounday constant as defined in Eqn #34 

* rs - Reference bounday radius of curvature as defined in Eqn #33 

* dmrdmrst - Derivative of mr wrt mrstar for mrpr calculation 

* dmrstdhr - Derivative of mrstar wrt hr for mrpr calculation 


* y ~ y location of point where dydx = 0. Defined in Eqn #llb, #17b 

* X - X location of point where dydx = 0. Defined in Eqn #lla, #17a 

* xihigh - The upper bounds for bisection method in xi 

* xilow - The lower bounds for bisection method in xi 

* ddcl - The constant used in the calculation of mrpr 

* gam - The ratio of specific heats 

* gmm - gam - 1. 

* gmp - gam + 1. 

* eps - The convergence criteria for the iteration. 

* 

* The first part of the subroutine sets all of the variables used 

* in this subroutine and prints that the program is seraching for 

* the eta that corresponds to the wall contour. 

* 


subroutine etawall 

real gam, cl, rs, xihigh, xilow 

real eps, xslO, mrstthroat, mstth 

external yminm,rtbis ! declares it an external function, needed for NR call 
logical succes, Isuper 

comroon/sonic/lsuper ! tells subroutine if mrst (Not mst!) in sonic area 
common/ input/ rs, cl, gam ! input values for passing to subroutines 
common/ throat/etaw, xithroat, mrstthroat , xslO 

* !save for later comparisons, value of eta 

* !at the nozzle wall, y=l. sps 

* remember that mrst depends only on xi, not on eta. So mrst = mrst 

* at same xi 

data tiny/2 . Oe-7/ 

data eps/l.Oe-5/ ! reduced from 0.00035 by sps with new routine 6-2-93 

* 

write (*,*) 'Search for Eta on Contour' 

* 

* This part of the program sets the denominator for the exponential 

* and sets the limits on the maximum and minimum values for xi 

* for the iteration. Since the value for dydx = 0 is near the axis, 

* the high value needs to be small. If it is not, the program will 

* not give accurate results. As the radius of curvature gets smaller, 

* the high value must get closer to zero. Since the only xi that is 

* used is the squared value, the low value should always be left at 

* zero to avoid division by zero in the dmrstdhr term when xi = 0. 

* 

xihigh = -0.03 
xilow = -0.2 

* 

* Call Numerical Recipes Routine to bracket for the root 

call zbrac (yminm, xilow, xihigh, succes) 
if (succes) then 

write (*,*) 'ETAWALL; successfully bracketed root in ', 

> xilow, xihigh 

else 

stop 'ETAWALL: failed to bracket root' 
end if 

* 

* Call Numerical Recipes Routine to iterate via Bisection for root 


xithl = rtbis (yminm, xilow, xihigh, eps) 

* xithl should be xi at the throat to within eps, now 

* function yminm also sets the value of parameters in common block /throat/ 

* BUT, last call from rtbis may not be at converged value. So make 

* this last call: 


ythroat = ymirnn(xithl) 

write (*,*) 'ETAWALL: ythroat converged to: ythroat 


* 

if (abs (xithl-xithroat) .gt. tiny) then 

write (*,*) 'ETAWALL: xithl , xithroat= ' ,xithl,xithroat 
stop 'ETAWALL: passing error' 
end if 

call getxy (mrstthroat , etaw, xithl , xth, yth,mstth) 
write (*,*) '***********************************' 

write (*,*) 'ETAWALL: set common block /throat/:' 
write (*,*) 'etaw=', etaw,' xithroat= ',xithroat 
write (*,*) 'mrstthroat= ', mrstthroat, ' xslO= ',xslO 
write (*,*) 'GETXY gives throat value for rostth= ',mstth 
write (*,*) 'and x,y at physical throat, wall= ',xth,yth 
write I*,*) '***********************************' 

return 
end 

* FUNCTION GETX 

* this is a function for ETAWALL that contains iterated function 

* part of the solver for the throat values of eta and xi. 

* modified from YMINM by sps 

* added sps 6-2-93 to allow use of Numerical Recipes bisector code 

* 

function getx(xi) 

real eterm, eden, hr, hrpr, hrdpr, etaw, ylpr, y3pr, yl 
real hrtpr, gmm, gmp, gam, cl, rs, y3, xi 
real xslO, mrst, mr, dmrdmrst, ddcl, dmrstdhr, mrpr 
real mrstthroat ,mrstpass 

logical Isuper Itrue if iterating for mrst (Not mstl) in supersonic region 
common/ sonic/ 1 super 

common/ input/rs, cl, gam ! input values for passing to subroutines 
common/throat/etaw,xithroat, mrstthroat, xslO lvalues at throat 
common/getxpass/xpass , etapass , mrstpass 

* (only used to pass to this subroutine) 

* 

gmm = gam-1.0 
gmp = gam+1.0 

* 

eden = 2.*rs*cl 

* 

* set sonic branch, see top rh column p. 1339 

if (xi .It. 0) then 

* Imrst subsonic, xi=0 is sonic line, exactly, for mrst (not mst!) 
Isuper = .false. 

else 

Isuper = .true, 
end if 

* Note that fig. 4, others, show supersonic flow at throat. 

* BUT THIS IS MST, NOT!! MRST WHICH IS WHAT IS BEING DECIDED HERE! 

* write (*,*) 'GETX: Isuper = ', Isuper, ', if this false, problem?' 

* 

eterm = EXP(-xi**2/eden) 
hr = cl*(l-eterm) + 1. 
hrpr = (xi/rs) *eterm 

hrdpr = (l./rs) *eterm - (xi**2/ (rs**2*cl) ) *eterm 

hrtpr = -(3.*xi/(rs**2*cl) ) *eterm + (xi**3 ./ (rs**3*cl**2 ) ) * 

+ eterm 

* 

* Call subroutine to determine the value of mrstar given hr. 


* 

call itermrst (hr , mrst) 

* Recursive call to rtbis inside itermrst! rename this call rtbis2 ! ! 

* 

* Calculate the value of mr and the needed derivatives. 

* 

mr = SQRT(2 . *mrst**2/ (gmp - gmm*mrst**2) ) 
dmrdmrst = (2 . *mrst*gmp/ (mr* (gmp - gmm*mrst**2) **2) ) 
ddcl = 0.5*gmp - 0 . 5*gmm*mrst**2 
if(mrst .LT. 0.99 .OR. mrst .GT. 1.01) then 
drarstdhr = 2 ./ (mrst*hr/ddcl - hr/mrst) 
mrpr = dmrdmrst*dmrstdhr*hrpr 
else 

mrpr = dmrdmrst*SQRT(2 ./ (gmp*rs) ) 
endif 

* 

* Calculate the values for the constants in the y equation (#llb, 

* #17b) and their derivatives. 

* 

yl = hr 

y3 = (hr/8 . 0) * (hr*hrdpr* (mr**2-l) - hrpr**2) 
ylpr = hrpr 

y3pr = (hrpr/8 . ) * (hr*hrdpr* (mr**2-l) - hrpr**2) + (hr/8.)* 

+ (hrpr*hrdpr* (mr**2-l) + hr*hrtpr* (mr**2-l) + 2*hr* 

+ hrdpr*mr*mrpr - 2*hrpr*hrdpr) 

* 

* 

* Now compute some other things also will be needed if this is root 

theta3 = (1./4 . ) * (hrdpr* (hr*hrpr*mr**2 + hr**2*mr*mrpr) ) + 

> (l./8.)*hr**2*hrtpr*(mr**2-l.) + (1./24 . ) *hrpr**3 
x2 = 0.5*hr*hrpr 

x4 = (3 ./32) * (hr**2*hrpr*hrdpr* (mr**2-l) ) - (1./96. ) * (hr* 

> hrpr**3) + theta3*hr*0 . 25 

X = xi - x2*etapass**2 - x4*etapass**4 
X = X - xslO ! shift so x=0 at throat 

getx = x-xpass ! return a value that should be zero, for rtbis 
mrstpass = mrst 

* write (*,5) xi , x, getx, mrstpass 

* 5 format (' GETX: return xi , x, getx, mrstpass= ^ , 4 (f 10 . 6 , lx) ) 

* write (*,*) 'and xpass was ',xpass 

* 

return 

end 


* FUNCTION YMINM 

* this is a function for ETAWALL that contains iterated function 

* part of the solver for the throat values of eta and xi. 

* modified from ETAWALL by sps 

* added sps 6-2-93 to allow use of Numerical Recipes bisector code 

* 

function yminm(xi) 

real eterm, eden, hr, hrpr, hrdpr, etaw, ylpr, y3pr, yl 
real hrtpr, gmm, gmp, gam, cl, rs, y3, xi 
real xslO, mrst, mr, dmrdmrst, ddcl, dmrstdhr, mrpr 
real mrstthroat 

logical Isuper Itrue if iterating for mrst in supersonic region 
common/ son ic/1 super 

common/input/rs,cl,gam ! input values for passing to subroutines 
common/throat/etaw,xi throat, mrstthroat, xslO lvalues at throat 



gmin = gam- 1.0 
gmp = gam+ 1 . 0 


eden = 2.*rs*cl 

★ 

* set sonic branch, see top rh column p. 1339, Note this is mrst, not mst!! 

if (xi .It. 0.0) then 1 subsonic, xi=0 is sonic line for msrt 
Isuper = .false, 
else 

Isuper = .true, 
end if 

* Note that fig. 4, others, show supersonic flow at throat. But this 

* is mst. Here we are looking at mrst, which is set by xi alone 

* write (*,*) 'YMINM: Isuper = ', Isuper ,', if this false, problem?' 

* 

eterm = EXP(-xi**2/eden) 
hr = cl*(l-eterm) + 1. 
hrpr = (xi/rs) *eterm 

hrdpr = (l./rs) *eterm - (xi**2/ (rs**2*cl) ) *eterm 

hrtpr = - ( 3 . *xi/ (rs**2*cl) ) *eterm + (xi**3 ./ (rs**3*cl**2) ) * 

+ eterm 

* 

* Call subroutine to determine the value of mrstar given hr. 

* 

call itermrst(hr, mrst) 

* Recursive call to rtbis inside itermrst! rename this call rtbis2!! 

* 

* Calculate the value of mr and the needed derivatives. 

* 

mr = SQRT(2 . *mrst**2/ (gmp - gmm*mrst**2) ) 
dmrdmrst = (2 . *mrst*gmp/ (mr* (gmp - gmm*mrst**2) **2) ) 
ddcl = 0.5*gmp - 0. 5*gmm*mrst**2 
if(mrst .LT. 0.99 .OR. mrst .GT. 1.01) then 
dmrstdhr = 2 ./ (mrst*hr/ddcl - hr/mrst) 
mrpr = dmrdmrst*dmrstdhr*hrpr 
else 

mrpr = dmrdmrst*SQRT(2 ./ (gmp*rs) ) 
endif 

* 

* Calculate the values for the constants in the y equation (#llb, 

* #17b) and their derivatives. 

* 

yl = hr 

y3 = (hr/8 . 0) * (hr*hrdpr* (mr**2-l) - hrpr**2) 
ylpr = hrpr 

y3pr = (hrpr/8 . ) * (hr*hrdpr* (mr**2-l) - hrpr**2) + (hr/8.)* 

+ (hrpr*hrdpr* (mr**2-l) + hr*hrtpr* (mr**2-l) + 2*hr* 

+ hrdpr*mr*mrpr - 2*hrpr*hrdpr) 

* 

* Set the value for eta where dydx = 0 as given by equation #25 and 

* calculate the value for y at that location. 

* 

argeta = -ylpr/y3pr 
if (argeta .gt. 0.0) then 
etaw = SQRT(-ylpr/y3pr) 
else 

write (*,*) 'YMINM; argeta = argeta,' It 0' 

write (*,*) 'for xi = ',xi 

stop 'fatal error, need sqrt of this' 


end if 

yminm = yl*etaw + y3*etaw**3 - 1.0 

* program converges if ymin =1.0, so subtract 1 here to provide zero 

* seeking function to Numerical Recipes 

* write out for debug 

* write (*,5) xi , etaw, yminm 

* 5 format (' YMINM: return xi , etaw, ymin-l= ' , 3 ( f 10 . 6 , lx) ) 

* 

* Now compute some other things also will be needed if this is root 

theta3 = (1./4 . ) * (hrdpr* (hr*hrpr*mr**2 + hr**2*mr*mrpr) ) + 

> (1 ./8 . ) *hr**2*hrtpr* (mr**2-l. ) + ( 1 ./24 . ) *hrpr**3 
x2 = 0.5*hr*hrpr 

x4 = (3 ./32) * (hr**2*hrpr*hrdpr* (mr**2-l) ) - (1./96. ) * (hr* 

> hrpr**3) + theta3*hr*0 . 25 

xslO = xi - x2*etaw**2 - x4*etaw**4 

xithroat = xi ! these two added sps for continuation later 
mrstthroat = mrst 

* 

return 

end 


* SUBROUTINE MRSTBAE 

* finds the mrst at the beginning angle and ending x 

subroutine mrstbae 

real mrstbeg , mrstend , mrstnow , mstout , mrstpass 
logical succes, Isuper 

external anglezer, rtbis, golden, mnbrak,getx 
common/sonic/ Isuper 

common/ input/rs, cl, gam ! input values for passing to subroutines 
common/ends/mr s tbeg , mrstend , thetaa , xend , thetanow , mrstnow , xbeg 
common/throat/etaw, xithroat, mrstthroat, xslO ! results at wall, throat 

* first iteration for continuation. Note that x,y = 0,1 at throat! 

* remember that mrst depends only on xi, not on eta! 

common/getxpass/xpass , etapass , mrstpass 

* (only used to pass to getx) 

data eps/l.Oe-5/ 

* data nplt/50/ 

data small/0. 01/! entry angle, closeness to specified value 

* 

thetanow = thetaa !do entry first 
xilow = -10.0 
xihigh = 0.0 

Isuper = .false. ! subsonic area, needed for mrst branch 
anglow = anglezer (xilow) 
anghigh = anglezer (xihigh) 
if (anglow*anghigh .It. 0.0) then 
succes = . true . 
else 

succes = .false, 
end if 

* remove zbrac, takes off too far looking for bracket 

* call zbrac(anglezer, xilow, xihigh, succes) 
if (succes) then 

write (*,*) 'MRSTBAE: bracketed entry angle, xi in ', xilow, xihigh 
write (*,*) 'corresponding angles are; ', 

> anglezer (xilow) +thetanow, anglezer (xihigh) +thetanow 

xibeg = rtbis (anglezer, xilow, xihigh, eps) 

* make one more call to set up final params 


from 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

* 


* 


aangzero = anglezer (xibeg) llast call to set up other params 
write (*,*) 'MRSTBAE: entry angle conv. within: aangzero 

mrstbeg = mrstnow 

write (*,*) 'MRSTBAE: xibeg, mrstbeg* xibeg, mrstbeg 
else I go for the minimum 

write (*,*) 'xilow,xihigh guessed as ' ,xilow, xihigh 
write (*,*) 'corresp. angles are: 

> anglezer (xilow) +thetanow, anglezer (xihigh) +thetanow 
write (*,*) 'entry angle is ', thetaa 

write (*,*) 'MRSTBAE: failed to bracket, find minimum' 

if (xilow .It. -100) xilow=-20 Idon't go way out 

delxi = (xihigh-xilow) /float (nplt-1) 

write (*,*) 'writing plot data to "entryang.deb” ' 

open (unit=9, file='entryang.deb') 

write (9,*) 'table of xi,angle(xi) : ' 

do 50 i = l,nplt 

xi = xilow+ ( i-1) *delxi 
write (9,*) xi , anglezer (xi) +thetanow 
50 continue 

close (unit=9) 

Call numerical recipes routine to find the minimum 
ax = -8.0 

bx = -5.0 !a guess near minimum 
cx = -3.0 

if (anglezer (bx) .It. anglezer (ax) .and. 

> anglezer (bx) .It. anglezer (cx) ) then 
write (*,*) 'minimum bracketed' 

else 

call mnbrak (ax, bx,cx, fa, fb, fc, anglezer) 
end if 

anglemin = golden (ax, bx,cx, anglezer, eps,ximin)+thetanow 
write (*,*) 'minimum entry angle is ', anglemin 
write (*,*) 'located at ximin* ',ximin 
if ( (anglemin-thetaa) /thetaa .It. small) then 

write (*,*) 'close to entry angle, keep this one' 
mrstbeg = mrstnow 
xibeg = ximin 

write (*,*) 'MRSTBAE: xibeg, mrstbeg* ', xibeg, mrstbeg 
else 

stop 'MRSTBAE: failed to find entry angle' 
end if 
end if 

call getxy ( mrstbeg , et aw , x iout , xout , yout , ms tout ) 
xbeg = xout 

write (*,*) 'MRSTBAE: xbeg* ',xbeg 

write (*,*) 'MRSTBAE: now finding exit xi on contour: ' 

etapass = etaw ! passing to getx 
xpass = xend 
xilow * -1.0 
xihigh = 8.0 
xzlow = getx (xilow) 
xzhigh = getx (xihigh) 
if (xzlow*xzhigh .It. 0.0) then 
succes = .true, 
else 

succes * .false, 
end if 

remove zbrac, goes off too far 


* call zbrac(getx,xilow,xihigh,succes) 
if (succes) then 

write (*,*) 'MRSTBAE: bracketed xend with xi in ' , xilow, xihigh 
else 

write (*,*) 'MRSTBAE: xend,xzhigh,xzlow= ' , xend, xzhigh, xzlow 
write (*,*) ' and xilow, xihigh= ' ,xilow,xihigh 

stop 'MRSTBAE; failed to bracket xend' 
end if 

xiend = rtbis(getx, xilow, xihigh,eps) 

* make one more call to set up final params 

xzero = getx (xiend) Hast call to set up other params 
write (*,*) 'MRSTBAE: xend conv. within: ', xzero 
mrstend = mrstpass ! passed back from getx 
write (*,*) 'MRSTBAE: xiend, mrstend= ', xiend, mrstend 

* 

return 

end 

* FUNCTION ANGLEZER 

* adapted by sps from MRSTBAE to contain function for Numerical Recipes 

* iteration, 6-2-93 

* finds the mrst at the beginning and ending angles 

function anglezer(xi) 

real thetaa, cl, rs, gam, gmp, gmm, etaw 
real mrstbeg, mrstend, mrstnow 

real xi, eterm, eden, hr, hrpr, hrdpr, thetacalc 
real hrtpr, mrst, mr, dmrdmrst, ddcl, dmrstdhr, rorpr 
real hrfpr, mdcl, mdc2, mdc3, mdc4, mdclpr, theta3 
real mrdpr, mdc2pr, mdc3pr, mdc4pr 

logical Isuper !true if iterating for mrst in supersonic region 
common/sonic/lsuper 

common/ input/rs, cl, gam ! input values for passing to subroutines 
common/ends/mrstbeg, mrstend, thetaa, xend, thetanow, mrstnow , xbeg 
common/throat/etaw, xithroat,mrstthroat,xslO Iresults at wall, throat from 

* first iteration for continuation. Note that x,y = 0,1 at throat! 

* remember that mrst depends only on xi, not on eta! 

* 

gmm = gam-1.0 

gmp = gam+1.0 

eden = 2.*rs*cl 

eterm = EXP(-xi**2/eden) 

hr = cl*(l-eterm) + 1. 

hrpr = (xi/rs) *eterm 

hrdpr = (l./rs) *eterm - (xi**2/ (rs**2*cl) ) *eterm 

hrtpr = -(3.*xi/(rs**2*cl))*eterm + (xi**3 ./ (rs**3*cl**2) ) * 

+ eterm 

hrfpr = -(3./(rs**2*cl) ) *eterm + (6.*xi**2/(rs**3*cl**2) ) * 

+ eterm - (xi**4/ (rs**4*cl**3) ) *eterm 

call itermrst(hr, mrst) ! finds mrst given hr 
mr = SQRT(2*mrst**2/ (gmp - gmm*mrst**2) ) 
ddcl = 0.5*gmp - 0. 5*gmm*mrst**2 
if (mrst .LT. 0.99 .OR. mrst .GT. 1.01) then 

dmrdmrst = (2 . *mrst*gmp/ (mr* (gmp - gmm*mrst**2) **2) ) 
dmrstdhr = 2 ./ (mrst*hr/ddcl - hr/mrst) 
mrpr = dmrdmrst*dmrstdhr*hrpr 
mrstpr = mrpr/ dmrdmrst 

ddmrdmrst = (4 . *mrstpr*gmp*ddcl - 4 . *mrst*gmp* (mrpr*ddcl - 
+ 2 . *mr*gmm*mrst*mrstpr) ) / (mr**2* (2 . *ddcl) **3) 

mdcl = gmp*mr 


+ + 


mdclpr = gTnp*mrpr 

mdc2= mrst**2 

mdc2pr = 2 . *mrst*mrstpr 

indc3 = hr 

mdc3pr = hrpr 

mdc4 = hr/mr**2 

mdc4pr = (hrpr*mr - 2 . *hr*mrpr) / (mr**3) 

mrdpr = ( (mdclpr* (mdc2* (mdc3-mdc4) ) -mdcl* (mdc2pr* (mdc3-mdc4) + 

+ mdc2* (mdc3pr-mdc4pr) ) )/ ( (mdc2* (mdc3-mdc4) ) **2) ) *hrpr+ 

+ dmrdmrst*dmrstdhr*hrdpr 

else 

mrpr = dmrdmrst*SQRT(2 ./ (<P^P*rs) ) 
mrstpr = mrpr/dmrdmrst 

ddmrdmrst = (4 . *mrstpr*gmp*ddcl - 4 . *mrst*gmp* (mrpr*ddcl - 
+ 2 . *mr*gmm*mrst*mrstpr) ) / (ror**2* (2 . *ddcl) **3) 

mrdpr = ddmrdmrst *SQRT (2 ./ (gmp*rs) ) 
endif 

theta3 = ( 1 ./ 4 • ) * (hrdpr* (hr*hrpr*mr**2 + hr**2*mr*mrpr) ) + 

+ ( 1 ./ 8 . ) *hr**2*hrtpr*(mr**2-l. ) + (1./24 . ) *hrpr**3 

thcl = hrdpr 
thclpr = hrtpr 
thc2 = hr*hrpr*mr**2 

thc2pr = (hrpr*mr)**2 + hr*hrdpr*mr**2 + 2*hr*hrpr*mr*mrpr 
thc3 = hr**2*mr*mrpr 

thc3pr = 2*hr*hrpr*mr*mrpr + (hr*mrpr) **2 + hr**2*mr*mrdpr 
thc4 = hr**2*hrtpr 

thc4pr = 2*hr*hrpr*hrtpr + hr**2*hrfpr 
thc5 = mr**2-l 
thcSpr = 2*mr*mrpr 
thc6pr = hrpr**2*hrdpr 

theta3pr = 0 . 25* (thclpr* (thc2+thc3) + thcl* (thc2pr+thc3pr) ) + 

+ 0. 125* (thc4pr*thc5 + thc5pr*thc4 + thc6pr) 

x2pr = (0. 5* (hrpr**2 + hr*hrdpr) ) 

x4cl = (3./32.)*((mr**2-l)*(2.*hr*hrpr**2*hrdpr +hr**2* 

+ hrdpr** 2 +hr** 2 *hrpr*hrtpr)+ 2 .*hr** 2 *hrpr*hrdpr*mr*mrpr) 

x4c2 = (1./96. ) *(hrpr**4 + 3 . *hr*hrpr**2*hrdpr) 
x4c3 = 0.25*(theta3pr*hr + theta3*hrpr) 
x4pr = (x4cl - x4c2 + x4c3) 
ylpr = hrpr 

y3pr = (hrpr/8 . ) * (hr*hrdpr* (mr**2-l) - hrpr**2) + (hr/8.)* 

(hrpr*hrdpr* (mr**2-l) + hr*hrtpr* (inr**2-l) + 2*hr* 
hrdpr *mr*mrpr - 2*hrpr*hrdpr) 
thetacalc = ATAN ( (ylpr*etaw+y3pr*etaw**3)/ (1. -x2pr*etaw**2 - 
+ x4pr*etaw**4) ) 

anglezer = thetacalc - thetanow !use thetanow for current sought value 
mrstnow = mrst 

write (*,*) 'ANGLEZER: return mrstnow, anglezer= 

> mrstnow , ang 1 e z e r 

write (*,*) 'and thetacalc, thetanow= ', thetacalc, thetanow 

return 

end 


SUBROUTINE ITERMRST 


The purpose of this subroutine is to find a value for mrstar given 
a value for hr. This is done by using a bisection method to 
iterate on mrstar to find the value of hr that is sent to the 



subroutine. 


The variables used in this subroutine are as follows; 

hr - The value of hr sent to this subroutine that is being matched. 

hrcalc - The calculated value of hr for a given inrstar. 

rorst - The value for mrstar that solves (13), sought here. 

gam - The ratio of specific heats. 

gmp - Gamma + 1 . 

gmm - Gamma - 1 . 

mrsthigh - The upper guess on mrst 
mrstlow - The lower guess on mrst 

The first part of the subroutine defines all of the variables and 
sets the guesses for mrst. 

subroutine itermrst(hr, mrst) 

real hrzero, hr, mrst, mrsthigh, mrstlow, mrstmax 
logical succes, Isuper Itrue if supersonic (2 branches!!) 
external hrzero, rtbis2 
comraon/sonic/1 super 

common /hriter/ hrwant ! a common block to pass the desired value 
common/ input/ rs, cl, gam ! input values for passing to subroutines 
data eps/l.Oe-6/ 

if hr = 1.0 exactly, handle as a special case to avoid singularities 
if (hr .eq. 1.0) then 
mrst = 1.0 

write (*,*) 'ITERMRST: returning exact sonic case' 
return 
end if 


gmm = gam-1.0 
gmp = gam+1.0 
mrstmax = sqrt (gmp/gmm) 

hrwant = hr 

if (Isuper) then !this is Isuper for mrst. Not local sonic (mst) 

mrstlow = 1.0 + eps 
mrsthigh = mrstmax - eps 
else 

mrstlow = eps 
mrsthigh =1.0- eps 
end if 

Call copies of numerical recipes routines to avoid recursive use, 
not supported in ordinary fortran 

call zbrac2 (hr zero, mrstlow, mrsthigh, succes) ! avoid, searches too far 

hrzlow = hrzero (mrstlow) 
hrzhigh = hrzero (mrsthigh) 
if (hrzlow*hrzhigh .It. 0.0) then 
succes = .true, 
else 

succes = .false, 
end if 

if (succes) then 

write (*,*) 'ITEMMRST: hr, Isuper = ', hr, Isuper 
write (*,*) 'ITERMRST; bracketed mrst in mrstlow, mrsthigh 
else 

write (*,*) 'mrstlow, mrsthigh= ', mrstlow, mrsthigh 
write (*,*) 'hrzlow, hrzhigh= ', hrzlow, hrzhigh 


write (*,*) 'hrwant= ',hrwant 
stop 'ITERMRST; failed to bracket root' 
end if 

itirst = rtbis 2 (hrzero,mrstlow,inrsthigh,eps) 

* write (*,*) 'ITERMRST: returned mrst= ',mrst,' at hr= ',hr, 

* > ' lsuper= ',lsuper 

* 

return 

end 

* following function added by sps to allow using Numerical Recipes routines 

function hrzero(mrst) 
real mrst 

common/hriter/hrwant !the desired value of hr 

common/ input/rs , cl , gam ! input values for passing to subroutines 

* 

gmro = gam-1.0 
gmp = gam+1.0 

argl = 0.5*gmp - 0 . 5*gmm*mrst**2 
if (argl .It. 0.0) then 

write (*,*) 'HRZERO: arg of exp is It 0' 
write (*,*) 'mrst,argl= mrst, argl 
stop 'fatal error' 
end if 

arg = 1 ./ (mrst* (argl) ** (1/gmm) ) 

if (arg .It. 0.0) stop 'HRZERO: arg of sqrt It zero' 

hrcalc = SQRT (arg) 

hrzero = hrcalc - hrwant 

return 

end 

* SUBROUTINE GETXY: 

* This subroutine modified from MSTAR 6-2-93 by sps. Given nozzle 

* definition, and guesses for mrst and eta, 

* it finds the x and y that correspond to the guesses. 

* Will be called 

* to find the mrst and eta for a particular 

* streamline needed for the square nozzle code 

* 1) this subroutine should only be called after the main program 

* has already been called, to initialize variables 

* 

* The purpose of this subroutine is to calculate Eqn #17 given 

* a specific value of mrstar on the axis and a value for a 

* streamline eta. 

* Then, at each point xi along the axis, 

* the values for hr and mr and their derivatives are calculated. 

* Using this information, eta is incremented from the axis (eta=0) 

* to the wall value (eta=etaw) . At each one of these points, the 

* corresponding value for x, y, theta, and mrstar is calculated. 

* 

* These are the variables used in this subroutine; 

* 

* hr - The value of h along the reference line (nozzle axis) . 

* Defined in Eqn. #12. 

* mr - The value of the Mach nuber along the reference line. 

* &&pr — The primed value of the function (either hr, mr, mrst, or theta3) 

* &&dpr - The second derivative of the function (either hr or mr) 

* hrtpr - The third derivative of hr 

* hrfpr - The fourth derivative of hr 


* xi - Velocity potential with units of length 

* etaw - The value of eta which corresponds to the wall 

* eterm - Epxonential of Eqn #22 

* x2 - Second term in series for x. Defined in Eqn #lla, #17a 

* x4 - Fourth term in series for x. Defined in Eqn #lla, #17a 

* yl ” First term in series for y. Defined in Eqn #llb, #17b 

* y3 ~ Third term in series for y. Defined in Eqn #llb, #17b 

* thetal - First term in series for theta. Defined in Eqn #llc, #17c 

* theta3 - Third term in series for theta. Defined in Eqn file, #17c 

* q2 - Second term in series for mstar. Defined in Eqn #lld, #17d 

* q4 - Fourth term in series for mstar. Defined in Eqn #lld, #17d 

* eta - The value of eta used from eta = 0 to eta = etaw. 

* mrst - The value of mstar on the reference line. 

* ra - The actual radius of curvature input by the user 

* cl - Reference bounday constant as defined in Eqn #34 

* rs - Reference bounday radius of curvature as defined in Eqn #33 

* thetaa - The actual initial flow angle of the wall input by 

* the user. 

* dmrdmrst - Derivative of mr wrt mrstar for mrpr calculation 

* dmrstdhr - Derivative of mrstar wrt hr for mrpr calculation 

* mrstbeg - The initial value of mrst determined from thetaa in mrstbae 

* dmrst - Delta mrst for the iteration in the direction of constant xi 

* eden - The denominator for the exponential te^ 

* deleta - Delta eta for the iteration in the direction of constant eta 

* ddmrdmrst - Derivative of dmrdmrst wrt xi for mrdpr calculation 

* ddmrstdhr - Derivative of dmrstdhr wrt xi for mrdpr calculation 

* xislO - The value of xi where dy/dx = 0 for shifting plots so minimum 

* point on top streamline is located at xi = 0. 

* ddcl - Constant for calculation of mr derivatives 

* ddc& - Constants 2-4 for calculation of ddmrstdhr 

* thc& - Constants 1-5 for calculation of theta3pr 

* imax - The maximum number of points in the xi direction 

* jmax - The maximum number of points in the eta direction 

* gam - The ratio of specific heats 

* gmm - gam - 1. 

* gmp - gam + 1 . 

* 

* Define the variables used in this subroutine, use common block 

* in order to find the eta which corresponds to the wall value. 

■k 

subroutine getxy (mrst , eta , xi , x , y , mst ) 

* given mrst and eta, returns xi, x, y, and mst 

real cl, rs, gmm, gmp, gam, mdc4pr 

real xi, hr, hrpr, hrdpr, hrtpr, hrfpr, q4, etaw 

real mr, dmrdmrst, mrst, eta, eterm, eden, mrstmax 

real ddcl, dmrstdhr, mrpr, mrstpr, ddmrdmrst, mdc2, mdc3 

real mdc4, mdcl, mdclpr, mdc2pr, mdc3pr, mrdpr, thcl, thc2, thc3 

real thc4, thc5, theta3pr, thetal, theta3, x2, x4, yl, y3, q2 

real x,y, theta, mst, mrstthroat 

real thclpr, thc2pr, thc3pr, thc4pr, theSpr, thc6pr 

common/throat/etaw, xithroat, mrstthroat, xslO !uses these, does not set 
common/ input/rs , cl , gam 

* 

gmm = gam-1.0 
gmp = gam+1.0 
mrstmax = sqrt (gmp/gmm) 

* 

if (etaw .eg. 0.0) stop ^GETSTR: etaw=0, fatal' 
eden = 2.*rs*cl 

* 



Calculate hr and its derivatives and mr and its derivatives. Also 
calculate the value of xi corresponding to the given value for 
mrstar from the do loop. 


* 

* 

* 

* 

* write (*,*) 'GETXY; received mrst,eta= ',mrst,eta 
if (mrst .ge. mrstmax) then 

write (*,*) 'GETXY: mrst,mrstmax= mrst, mrstmax 
stop 'GETXY: mrst ge mrstmax' 
end if 

if (mrst .le. 0.0) stop 'GETXY: mrst le zero' 

hr = SQRT( (l./mrst) * (0. 5*gmp - 0. 5*gmm*mrst**2) ** (-1 ./ (gmm) ) ) 
argxil = cl/ (cl - hr + 1.) 
if (argxil .ge. 0.0) then 

argxi = 2 . *rs*cl*LOG (argxil) 
else 1 singular 

write (*,*) 'GETXY: argxil= ', argxil,' mrst= ',mrst 
write (*,*) 'and hr= ',hr,' rs,cl= ',rs,cl 
stop 'this causes fatal error, argxil (of log) It 0.0' 
end if 

if (argxi .ge. 0.0) then 
xi = SQRT (argxi) 
else 

write (*,*) 'GETXY: argxi= ', argxi,' mrst= ',mrst 
write (*,*) 'and hr= ',hr,' rs,cl= ',rs,cl 
stop 'this causes fatal error, argxi It 0.0' 
end if 

if (mrst .LE. 1.) then 
xi = -xi 
endif 

eterm = EXP(-xi**2/eden) 
hrpr = (xi/rs) *eterm 

hrdpr = (1 ./rs) *eterm - (xi**2/ (rs**2*cl) ) *eterm 

hrtpr = - ( 3 . *xi/ (rs**2*cl) ) *eterm + (xi**3 ./ (rs**3*cl**2) ) * 

+ eterm 

hrfpr = - (3 ./ (rs**2*cl) ) *eterm + (6. *xi**2/ (rs**3*cl**2) ) * 

+ eterm - (xi**4/ (rs**4*cl**3) ) *eterm 

mr = SQRT(2*mrst**2/ (gmp - gmm*mrst**2) ) 
ddcl = 0.5*gmp - 0 . 5*gmm*mrst**2 

dmrdmrst = (2 . *mrst*gmp/ (mr* (gmp - gmm*mrst**2) **2) ) 
if (mrst .LT. 0.99 .OR. mrst .GT. 1.01) then 

* use special formula derived near singularity in derivatives 

dmrstdhr = 2 ./ (mrst*hr/ddcl - hr/mrst) 
mrpr = dmrdmrst*dmrstdhr*hrpr 
mrstpr = mrpr/ dmrdmrst 

ddmrdmrst = (4 . *mrstpr*gmp*ddcl - 4 . *mrst*gmp* (mrpr*ddcl - 
+ 2 . *mr*gmm*mrst*mrstpr) ) / (mr**2* (2 . *ddcl) **3) 

mdcl = gmp*mr 
mdclpr = gmp*mrpr 
mdc2= mrst**2 
mdc2pr = 2 . *mrst*mrstpr 
mdc3 = hr 
mdc3pr = hrpr 
mdc4 = hr/mr**2 

mdc4pr = (hrpr*mr - 2 . *hr*mrpr)/ (mr**3) 

mrdpr = ( (mdclpr* (mdc2* (mdc3-mdc4) ) -mdcl* (mdc2pr* (mdc3-mdc4) + 
+ mdc2* (mdc3pr-mdc4pr) ) )/ ( (mdc2* (mdc3-mdc4) ) **2) ) *hrpr+ 

+ dmrdmrst *dmrstdhr*hrdpr 

else 

mrpr = dmrdmrst*SQRT(2 ./ (gmp*rs) ) 
mrstpr = mrpr/dmrdmrst 


ddmrdmrst = (4 . *mrstpr*gmp*ddcl - 4 . *inrst*ginp* (mrpr*ddcl - 
+ 2 . *mr*gmm*mrst*inrstpr) ) / (inr**2* (2 . *ddcl) **3) 

mrdpr = ddmrdmrst*SQRT(2 ./ ) 
end if 

Calculate the values for thetal, theta3, theta3pr, yl, y3, x2, 
x4, q2, q4 from Eqn #17a-#17d. All of these remain constant at 
any eta for a given value of xi since they are only dependant 
on hr, mr, and their derivatives. 

thcl = hrdpr 
thclpr = hrtpr 
thc2 = hr*hrpr*mr**2 

thc2pr = (hrpr*mr) **2 + hr*hrdpr*mr**2 + 2*hr*hrpr*mr*mrpr 
thc3 = hr**2*mr*mrpr 

thc3pr = 2*hr*hrpr*mr*mrpr + (hr*mrpr)**2 + hr**2*mr*mrdpr 
thc4 = hr**2*hrtpr 

thc4pr = 2*hr*hrpr*hrtpr + hr**2*hrfpr 
thc5 = mr**2-l 
thc5pr = 2*mr*mrpr 
thc6pr = hrpr**2*hrdpr 

theta3pr = 0. 25* (thclpr* (thc2+thc3) + thcl* (thc2pr+thc3pr) ) + 

+ 0. 125* (thc4pr*thc5 + thc5pr*thc4 + thc6pr) 

thetal = hrpr 

theta3 = (1./4 .)* (hrdpr* (hr*hrpr*mr**2 + hr**2*mr*mrpr) ) + 

+ (1./8 . ) *hr**2*hrtpr* (mr**2-l. ) + (1./24 . ) *hrpr**3 

x2 = 0.5*hr*hrpr 

x4 = (3 ./32) * (hr**2*hrpr*hrdpr* (mr**2-l) ) - (1./96. ) * (hr* 

+ hrpr**3) + theta3*hr*0 . 25 

yl = hr 

y3 = (hr/8) * (hr*hrdpr* (mr**2-l) - hrpr**2) 
q2 = 0.5*hr*hrdpr 

q4 = 0.25*(hr**2*hrdpr**2) + (3 ./32 . ) * (hr**2*hrdpr**2* (mr**2- 
+ 1) ) + (1/32 . ) * (hr*hrdpr*hrpr**2) + 0. 25*hr*theta3pr 

x2pr = (0. 5* (hrpr**2 + hr*hrdpr) ) 

x4cl = ( 3 ./ 32 . ) * ( (mr**2-l) * (2 . *hr*hrpr**2*hrdpr +hr**2* 

+ hrdpr**2+hr**2*hrpr*hrtpr) +2 . *hr**2*hrpr*hrdpr*mr*mrpr) 

x4c2 = (1./96. ) * (hrpr**4 + 3 . *hr*hrpr**2*hrdpr) 
x4c3 = 0.25*(theta3pr*hr + theta3*hrpr) 
x4pr = (x4cl - X4c2 + X4c3) 
ylpr = hrpr 

y3pr = (hrpr/8 . ) * (hr*hrdpr* (mr**2-l) - hrpr**2) + (hr/8.)* 

+ (hrpr*hrdpr* (mr**2-l) + hr*hrtpr* (mr**2-l) + 2*hr* 

+ hrdpr*mr*mrpr - 2*hrpr*hrdpr) 

Calculate Eqn #17. 

Calculate Eqn #17a. Then, if mrstar is less than one, make the value 
for X the negative value for plotting purposes. This is because 
the value for xi is always positive along the axis. However, since 
the paper makes the assumption that the sonic line on the axis occurs 
at X = 0 (pp. 1339 paragraph before section #3, and figure #3), it 
is known that for all values of mrstar less than 1, the x value 
should be negative. 

X = xi - x2*eta**2 - x4*eta**4 

Shift the X axis by xslO so that the point where x = 0 corresponds 
to the nozzle throat. 



* 

* 

* 


X 


X - xslo 


Calculate Eqn #17b. 
y = yl*eta + y3*eta**3 

* 

* Calculate Eqn #17c. 

* 

theta = ATAN( (ylpr*eta+y3pr*eta**3)/(l.-x2pr*eta**2 - 
+ x4pr*eta**4) ) 

* 

* Calculate Eqn #17d. 

* 

* this gives the local mach number, which depends on eta as well 

* as on xi. 

mst = mrst*(l + q2*eta**2 + q4*eta**4) 

* 

return 

end 



